TH1: Geospatial Analytics for Social Good

Take Home Exercise 1: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar

Author

Kendrick Teo

Published

September 7, 2024

Modified

September 21, 2024

TH1.1 Setting the Scene

Millions of people have their lives shattered by armed conflict every year. One of these is the Myanmar Civil War, a significant escalation of the long-running Myanmar Conflict in response to the 2021 coup d’etat.

Geospatial analtyics holds the potential to address complex problems facing society, such as this one. This study serves to discover the sptial and spatio-temporal distribution (spread) of the latest armed conflict in Myanmar by applying spatial point pattern analysis (SPPA) methods.

TH1.1.1 Loading R packages

The R packages we will use today are:

  • sf
  • tmap
  • spatstat
  • sparr
  • raster
  • maptools
pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse)

TH1.2 The Data

Armed conflict data in Myanmar between 2010 and 2023 was downloaded from the Armed Conflict Location & Event Database (ACLED), an independent, impartial, international non-profit organization collecting data on violent conflicts and protests in all countries and territories in the world (Raleigh et al., 2023). We will be superimposing these locations with the geogrpahical boundary and subdivisions of the country, from the Myanmar Information Management Unit (MIMU).

We are interested in conflict events from 2019 onwards - starting from a year before the COVID-19 pandemic and 2 years before the 2021 coup - cumulating into the civil war of today.

TH1.2.1 Loading armed conflict data into tibble object

Here we load our aspatial armed conflict data into a tibble object.

acled_mya <- read_csv("data/2021-01-01-2024-09-16-Southeast_Asia-Myanmar.csv")
Rows: 53653 Columns: 31
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (20): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (11): year, time_precision, inter1, inter2, interaction, iso, latitude, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(acled_mya)
# A tibble: 6 × 31
  event_id_cnty event_date         year time_precision disorder_type  event_type
  <chr>         <chr>             <dbl>          <dbl> <chr>          <chr>     
1 MMR66329      06 September 2024  2024              1 Political vio… Battles   
2 MMR66334      06 September 2024  2024              1 Political vio… Explosion…
3 MMR66353      06 September 2024  2024              1 Political vio… Explosion…
4 MMR66365      06 September 2024  2024              1 Political vio… Explosion…
5 MMR66368      06 September 2024  2024              1 Political vio… Battles   
6 MMR66422      06 September 2024  2024              1 Political vio… Explosion…
# ℹ 25 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
#   inter1 <dbl>, actor2 <chr>, assoc_actor_2 <chr>, inter2 <dbl>,
#   interaction <dbl>, civilian_targeting <chr>, iso <dbl>, region <chr>,
#   country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>, location <chr>,
#   latitude <dbl>, longitude <dbl>, geo_precision <dbl>, source <chr>,
#   source_scale <chr>, notes <chr>, fatalities <dbl>, tags <chr>,
#   timestamp <dbl>

TH1.2.2 Loading administrative boundaries into sf object

Three representations of Myanmar’s geography exist on the MIMU repository - Admin1 subdivides the country into its states and regions only, while Admin2 subdivides the country by its smaller districts, and Admin3 its townships. Further, the ACLED labels each incident with all 3 representations. For simplicity and ease of understanding, we shall use the Admin1 representation.

mmr_admin1 <- st_read(dsn = "data/mmr_polbnda_adm1_250k_mimu_1", layer = "mmr_polbnda_adm1_250k_mimu_1")
Reading layer `mmr_polbnda_adm1_250k_mimu_1' from data source 
  `/Users/kendricktty/Gits/smu_cs/is415-site/TakeHome/TakeHome1/data/mmr_polbnda_adm1_250k_mimu_1' 
  using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
mmr_admin1
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4
                         geometry
1  MULTIPOLYGON (((95.20798 15...
2  MULTIPOLYGON (((96.17964 19...
3  MULTIPOLYGON (((93.36931 24...
4  MULTIPOLYGON (((97.59674 28...
5  MULTIPOLYGON (((97.1759 19....
6  MULTIPOLYGON (((97.81508 16...
7  MULTIPOLYGON (((94.11699 22...
8  MULTIPOLYGON (((96.14023 23...
9  MULTIPOLYGON (((97.73689 15...
10 MULTIPOLYGON (((96.32013 20...
qtm(mmr_admin1) # Plots the national boundaries of Myanmar

TH1.2.3 Creating sf data frame from aspatial data

The ACLED tibble contains coordinates, making it useful for plotting on our map as points. We can therefore use it to create an sf data frame using which we can plot our points on a map. The EPSG format of the import coordinates should be 4326, corresponding to the WGS84 Geographic Coordinate System.

acled_mya_sf <- st_as_sf(acled_mya, coords = c("longitude", "latitude"), crs = 4326)

TH1.3 A brief summary of Myanmar and the Myanmar Conflict

The second largest country in Southeast Asia at 676,578 square kilometres and an ASEAN member state, Myanmar features fertile tropical deltas in the south and the rugged Himalayan foothills to its north. It shares borders with China to the north and northeast, Laos and Thailand to the east and Southeast, and Bangladesh and India to the west and northwest. Despite having a 2800km coastline providing access to deep-sea ports and an abundance in natural resources including arable land, forests, minerals and natural gas, conflict-ridden Myanmar is considered a least developed country, with a GNP per capita of US$1144 (2011) (Myanmar Information Management Unit, n.d.) and a HDI of 0.6 (UNDP, 2024). Consequently, agriculture accounts for 36% of GDP and 60-70% of employment.

Myanmar is organised into seven states, seven regions, and one union territory - its capital, Naypyitaw, located at around the country’s geographic centre. The states - Chin, Kachin, Kayah, Kyain, Mon, Rakhine and Shan - are largely populated by the national ethnic communities, while most of the inhabitants of its regions - Ayeyarwady, Bago, Mandalay, Sagaing, Tanintharyi and largest city Yangon - are of the majority Bamar (Burmese) ethnicity. The following code chunk plots the states and regions of Myanmar into an interactive map.

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(mmr_admin1) + tm_polygons()
tmap_mode('plot')
tmap mode set to plotting

Crucially, Myanmar is one of the world’s most ethnically diverse countries, with as many as 135 different ethnic groups boasting a rich tapestry of culture and religious history (Maizland, 2022). While this may paint a rosy picture of harmony on the surface, this is far from the case in the country. While the majority Bamar enjoy a privileged position in society and have held of government and military positions, many of the country’s minority groups grapple with systemic discrimination, lack of economic opportunity and development, minimum representation in government, and abuses at the hands of the military.

The code chunk below visualises Myanmar’s ethnic breakdown in a pie chart, based on data from the US CIA World Factbook and presented by Maizland (2022).

# Plot the composition of ethnic groups in Myanmar
x <- c(2, 2, 3, 4, 5, 7, 9, 68)
y <- c("Indian", "Mon", "Chinese", "Rakhine", "Other", "Karen", "Shan", "Bamar")

y <- paste(y, x)
y <- paste(y, "%", sep="")

pie(x, labels = y, main = "Composition of Ethnic Groups in Myanmar (Maizland, 2022)", col = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))

# legend("topright", y, cex = 0.8, fill = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"))

Because of this disparity and partially as a result of the divide-and-conquer strategy used by the British during colonial rule, lengthy armed conflicts have ensued since Myanmar’s independence between the ethnic groups (comprising ethnic armed organisations and smaller militias) and the Tatmadaw, the military, producing what has become known as the world’s longest continuing civil war (Rieffel, 2019), over issues from greater autonomy to control of natural resources. This did not stop even after the country underwent a democratic reform in the early 2010s - since 2017, the Tatmadaw had been mounting a brutal campaign against the ethnic Rohingya (Albert & Maizland, 2020). While the administration of then-state chancellor Aung San Suu Kyi defended their actions, restricting press freedoms in the process, rights groups and UN officials claim the Tatmadaw’s actions amounted to genocide.

TH1.4 Data wrangling

TH1.4.1 Coordinate system conversion

The EPSG area code for Myanmar is 4239, while the WGS 84-compatible code is 32647. Since we’re creating an owin object later, we need to first convert our coordinate system.

acled_mya_sf <- acled_mya_sf %>% st_transform(crs = 32647)
mmr_admin1 <- mmr_admin1 %>% st_transform(crs = 32647)
tm_shape(mmr_admin1) + tm_polygons() + tm_shape(acled_mya_sf) + tm_dots(size = 0.05)

TH1.4.2 Feature creation

Since the data is timestamped, we are able to plot and compare the frequency of conflict events over time. Before we do, though, we need to standardise the date stamps, then compartmentalise them into their year, month and day components, so that we can perform spatial-temporal point patterns analysis (STPPA) later. We will also aim to segregate events by whether or not they occurred before the coup on February 1, 2021.

A visualisation of the spatial-temporial distribution of points can be found in Annex A. Notice the consistent appearance throughout the entire study period of a large cluster of points towards the west of the country - roughly the location of Rahkine state on the border with Bangladesh.

Additionally, the ACLED dataset splits the conflict incidents in Shan state into North, South and East, and incidents in the Bago region into East and West. We need to combine these into their representations in the polygon set.

acled_mya_sf <- acled_mya_sf %>%
    mutate(event_date = dmy(event_date)) %>%
    mutate(DayOfYear = yday(event_date)) %>%
    mutate(Month_num = format(event_date, "%Y-%m")) %>%  # Year-Month format
    mutate(Month_fac = factor(format(event_date, "%B %Y"))) %>% # Full month
    mutate(Quarter_num = ((year(event_date) - 2021) * 4) + quarter(event_date)) %>% # Numbers quarters for STKDE
    mutate(Quarter = paste(year(event_date), "-", quarter(event_date), "Q"))  %>% # Quarter format
    mutate(admin1 = case_when(
      admin1 %in% c("Shan-North", "Shan-South", "Shan-East") ~ "Shan",
      admin1 %in% c("Bago-East", "Bago-West") ~ "Bago",
      TRUE ~ admin1
    ))
 [1] "Kachin"      "Shan"        "Magway"      "Sagaing"     "Chin"       
 [6] "Rakhine"     "Mandalay"    "Tanintharyi" "Kayah"       "Yangon"     
[11] "Kayin"       "Mon"         "Ayeyarwady"  "Bago"        "Nay Pyi Taw"
 [1] 15 14 13 12 11 10  9  8  7  6  5  4  3  2  1
[1] 15
 [1] "2024 - 3 Q" "2024 - 2 Q" "2024 - 1 Q" "2023 - 4 Q" "2023 - 3 Q"
 [6] "2023 - 2 Q" "2023 - 1 Q" "2022 - 4 Q" "2022 - 3 Q" "2022 - 2 Q"
[11] "2022 - 1 Q" "2021 - 4 Q" "2021 - 3 Q" "2021 - 2 Q" "2021 - 1 Q"
# Plot points by quarter
tm_shape(mmr_admin1) + tm_polygons() + tm_shape(acled_mya_sf) + tm_dots(size = 0.05) + tm_facets(by="Quarter_num", free.coords=FALSE, drop.units = TRUE)

Simple feature collection with 6 features and 34 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 297846.6 ymin: 1920825 xmax: 389155.2 ymax: 2653913
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 6 × 35
  event_id_cnty event_date  year time_precision disorder_type      event_type   
  <chr>         <date>     <dbl>          <dbl> <chr>              <chr>        
1 MMR11022      2021-01-02  2021              1 Political violence Explosions/R…
2 MMR10936      2021-01-02  2021              1 Political violence Battles      
3 MMR10933      2021-01-01  2021              1 Political violence Battles      
4 MMR10920      2021-01-01  2021              1 Political violence Violence aga…
5 MMR10906      2021-01-01  2021              1 Political violence Battles      
6 MMR10932      2021-01-01  2021              1 Political violence Battles      
# ℹ 29 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
#   inter1 <dbl>, actor2 <chr>, assoc_actor_2 <chr>, inter2 <dbl>,
#   interaction <dbl>, civilian_targeting <chr>, iso <dbl>, region <chr>,
#   country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>, location <chr>,
#   geo_precision <dbl>, source <chr>, source_scale <chr>, notes <chr>,
#   fatalities <dbl>, tags <chr>, timestamp <dbl>, geometry <POINT [m]>,
#   DayOfYear <dbl>, Month_num <chr>, Month_fac <fct>, Quarter_num <dbl>, …

TH1.4.3 Using owin to confine study area within Myanmar borders

Finally, we need to create an owin object that will be combined with any ppp objects we create later for SPPA.

mmr_admin1_owin <- as.owin(mmr_admin1)
plot(mmr_admin1_owin)

# summary(mmr_admin1_owin)

TH1.5 Preliminary analysis

TH1.5.1 Exploratory Data Analysis (EDA)

We can quickly gain some meaningful insights from our dataset using basic EDA techniques. For instance, with the following code chunk, we can retrieve the total number of fatalities from the armed violence during the study period on a quarterly basis, and print its corresponding summary statistics.

TH1.5.1.1 Getting quarterly distribution of fatalities and summary statistics

acled_mya_fatalities_by_quarter <- acled_mya_sf %>%
  group_by(Quarter_num) %>%
  summarise(`All fatalities` = sum(fatalities)) %>%
  ungroup()
# acled_mya_fatalities_by_quarter
summary(acled_mya_fatalities_by_quarter)
  Quarter_num   All fatalities          geometry 
 Min.   : 1.0   Min.   : 803   MULTIPOINT   :15  
 1st Qu.: 4.5   1st Qu.:3640   epsg:32647   : 0  
 Median : 8.0   Median :3955   +proj=utm ...: 0  
 Mean   : 8.0   Mean   :3935                     
 3rd Qu.:11.5   3rd Qu.:4798                     
 Max.   :15.0   Max.   :5862

We can establish that 3935 people die from the armed violence each quarter.

TH1.5.1.2 Plotting quarterly distribution of fatalities as a line chart

ggplot(acled_mya_fatalities_by_quarter, aes(x=Quarter_num, y=`All fatalities`)) +     
  # geom_bar(stat = "identity", fill = "red4") +
  geom_line(color = "red4") + geom_point() +
    labs(title = "Total Fatalities from the Armed Conflict in Myanmar by Quarter", subtitle = "1 = Q1 2021, 6 = Q2 2022, 11 = Q3 2023, 15 = Q3 2024",
         x = "Quarter",
         y = "All Fatalities") +
    theme_minimal()

There were fewer deaths from armed conflict in the first three quarters of 2021 - at the very start of the conflict - than Q4 2021 onwards.

TH1.5.1.3 Getting quarterly distribution of incidents involving civilian targeting

Using the same code structure, we can also plot the distribution of armed conflict events involving civilian targeting.

acled_mya_civilianTargeting_by_quarter <- acled_mya_sf %>%
  group_by(Quarter_num) %>%
  summarise(
    Count = n(),
    Civilians_targeted = sum(civilian_targeting == "Civilian targeting", na.rm = TRUE)
  ) %>%
  mutate(Percentage_incidents_civilians_targeting = (Civilians_targeted / Count) * 100) %>%
  ungroup()
# acled_mya_civilianTargeting_by_quarter
summary(acled_mya_civilianTargeting_by_quarter)
  Quarter_num       Count      Civilians_targeted          geometry 
 Min.   : 1.0   Min.   :2100   Min.   :354.0      MULTIPOINT   :15  
 1st Qu.: 4.5   1st Qu.:3229   1st Qu.:641.5      epsg:32647   : 0  
 Median : 8.0   Median :3689   Median :697.0      +proj=utm ...: 0  
 Mean   : 8.0   Mean   :3577   Mean   :681.8                        
 3rd Qu.:11.5   3rd Qu.:3999   3rd Qu.:758.0                        
 Max.   :15.0   Max.   :4727   Max.   :841.0                        
 Percentage_incidents_civilians_targeting
 Min.   : 9.581                          
 1st Qu.:18.388                          
 Median :20.279                          
 Mean   :19.629                          
 3rd Qu.:22.451                          
 Max.   :26.961                          
ggplot(acled_mya_civilianTargeting_by_quarter, aes(x=Quarter_num, y=Percentage_incidents_civilians_targeting)) +     
  # geom_bar(stat = "identity", fill = "red4") +
  geom_line(color = "skyblue1") + geom_point() +
    labs(title = "Percentage of Incidents involving civilian targeting in Myanmar by Quarter", subtitle = "1 = Q1 2021, 6 = Q2 2022, 11 = Q3 2023, 15 = Q3 2024",
         x = "Quarter",
         y = "% Incidents with civilian targeting") +
    theme_minimal()

The percentage of violent conflict incidents targeting civilians trended upwards throughout the study duration. On average, 19.6% of all armed conflict incidents throughout the study period targeted civilians.

TH1.5.1.4 Plotting proportion of fatalities by type of conflict

We can also plot the proportion of fatalities by the type of conflict. The ACLED classifies armed conflict events by mode and intensity of attack and by whether civilians are targeted. From the graph and pie chart below, the largest proportion of fatalities comes from battles, closely followed by explosions and remote violence. However, violence against civilians makes up a sizeable 11.4% of fatalities over the study period.

unique(acled_mya_sf$event_type)
[1] "Battles"                    "Explosions/Remote violence"
[3] "Strategic developments"     "Protests"                  
[5] "Violence against civilians" "Riots"                     
acled_mya_fatalities_by_event_type <- acled_mya_sf %>%
  group_by(event_type) %>%
  summarise(`All fatalities` = sum(fatalities)) %>%
  ungroup()
summary(acled_mya_fatalities_by_event_type)
  event_type        All fatalities             geometry
 Length:6           Min.   :  127.0   MULTIPOINT   :6  
 Class :character   1st Qu.:  288.8   epsg:32647   :0  
 Mode  :character   Median : 3556.0   +proj=utm ...:0  
                    Mean   : 9838.2                    
                    3rd Qu.:11037.5                    
                    Max.   :39059.0                    
ggplot(acled_mya_fatalities_by_event_type, aes(x=event_type, y=`All fatalities`)) +
  geom_bar(stat = "identity", fill = "red4") +
      labs(title = "Fatalities by event type in Myanmar",
         x = "Quarter",
         y = "All Fatalities") +
    theme_minimal()

total_fatalities <- sum(acled_mya_fatalities_by_event_type$`All fatalities`)
percentages <- 100 * acled_mya_fatalities_by_event_type$`All fatalities` / total_fatalities

# Set plotting area size
par(mar = c(1, 1, 2, 1))  # Adjust margins if needed

pie(acled_mya_fatalities_by_event_type$`All fatalities`, labels = paste0(round(percentages, 1), "%"), 
    col = rainbow(length(acled_mya_fatalities_by_event_type$event_type)), 
    main = "Fatalities by event type in Myanmar",
    cex = 1  
)
    
legend("topright", paste(acled_mya_fatalities_by_event_type$event_type, ":", round(percentages, 1), "%"), fill = rainbow(length(acled_mya_fatalities_by_event_type$event_type)), cex = 0.7)

We can further narrow down the scope of this analysis by limiting our data to a certain year, and to events where the civilian_targeting column is filled. The year 2022 saw the highest fatality numbers across the entire study period, so we will take a closer look.

acled_mya_fatalities_by_event_type_2022 <- acled_mya_sf %>%
  filter(year == "2022") %>%
  filter(civilian_targeting == "Civilian targeting") %>%
  group_by(event_type) %>%
  summarise(`All fatalities` = sum(fatalities)) %>%
  ungroup()

ggplot(acled_mya_fatalities_by_event_type_2022, aes(x=event_type, y=`All fatalities`)) +
  geom_bar(stat = "identity", fill = "red4") +
      labs(title = "Fatalities by event type in 2022 with civilian_targeting filled",
         x = "Quarter",
         y = "All Fatalities") +
    theme_minimal()

TH1.5.1.5 Mapping fatalities and violence against civilians

Finally, we can plot choropleth maps displaying 2 crucial pieces of information: the locations with the greatest number of fatalities, and the locations where civilians are targeted in the violence the most in 2022. In the code chunks we will write to create the plots, we have to refer to the original tibble objects instead of the sf objects.

acled_mya_fatalities_by_admin1_2022 <- acled_mya %>%
  filter(year == 2022) %>%
  mutate(admin1 = case_when(
      admin1 %in% c("Shan-North", "Shan-South", "Shan-East") ~ "Shan",
      admin1 %in% c("Bago-East", "Bago-West") ~ "Bago",
      TRUE ~ admin1
    )) %>%
  group_by(admin1) %>%
  summarise(`All_Fatalities` = sum(fatalities)) %>%
  ungroup()
combined <- left_join(mmr_admin1, acled_mya_fatalities_by_admin1_2022, by = c("ST" = "admin1"))
tm_shape(combined) + 
  tm_fill("All_Fatalities", n = 8, style = "kmeans") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(title = "admin1 by number of fatalities, 2022", legend.outside = TRUE)

The above map shows Sagaing region as the location with the highest number of fatalities in 2022, followed by Magaway region.

#|echo: false
combined
Simple feature collection with 15 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -210008.6 ymin: 1072026 xmax: 724647.6 ymax: 3158467
Projected CRS: WGS 84 / UTM zone 47N
First 10 features:
   OBJECTID          ST ST_PCODE           ST_RG          ST_MMR PCode_V
1         1  Ayeyarwady   MMR017          Region  ဧရာဝတီတိုင်းဒေသကြီး     9.4
2         2        Bago   MMR111          Region    ပဲခူးတိုင်းဒေသကြီး     9.4
3         4        Chin   MMR004           State       ချင်းပြည်နယ်     9.4
4         5      Kachin   MMR001           State       ကချင်ပြည်နယ်     9.4
5         6       Kayah   MMR002           State       ကယားပြည်နယ်     9.4
6         7       Kayin   MMR003           State        ကရင်ပြည်နယ်     9.4
7         8      Magway   MMR009          Region   မကွေးတိုင်းဒေသကြီး     9.4
8         9    Mandalay   MMR010          Region မန္တလေးတိုင်းဒေသကြီး     9.4
9        10         Mon   MMR011           State         မွန်ပြည်နယ်     9.4
10       11 Nay Pyi Taw   MMR018 Union Territory        နေပြည်တော်     9.4
   All_Fatalities                       geometry
1              57 MULTIPOLYGON (((93411.72 17...
2             364 MULTIPOLYGON (((203949.9 21...
3            1158 MULTIPOLYGON (((-72918.03 2...
4             756 MULTIPOLYGON (((362696.3 31...
5            1008 MULTIPOLYGON (((309155.7 22...
6            1021 MULTIPOLYGON (((373551.3 18...
7            2488 MULTIPOLYGON (((-1717.607 2...
8            1068 MULTIPOLYGON (((208184.3 26...
9             711 MULTIPOLYGON (((364238.4 16...
10             13 MULTIPOLYGON (((220155 2248...
acled_mya_civiliansTargeted_by_admin1_2022 <- acled_mya %>%
  filter(year == 2022) %>%
  mutate(admin1 = case_when(
      admin1 %in% c("Shan-North", "Shan-South", "Shan-East") ~ "Shan",
      admin1 %in% c("Bago-East", "Bago-West") ~ "Bago",
      TRUE ~ admin1
    )) %>%
  group_by(admin1) %>%
  summarise(
    count = n(),
    Civilians_targeted = sum(civilian_targeting == "Civilian targeting", na.rm = TRUE)
  )%>%
  ungroup() %>%
  mutate(`%_with_Civilian_Targeting` = (`Civilians_targeted` / `count`) * 100)
# acled_mya_civiliansTargeted_by_admin1_2022
combined <- left_join(mmr_admin1, acled_mya_civiliansTargeted_by_admin1_2022, by = c("ST" = "admin1"))
# combined
tm_shape(combined) + 
  tm_fill("%_with_Civilian_Targeting", n = 8, style = "kmeans") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(title = "admin1 by % of incidents targeting civilians, 2022", legend.outside = TRUE)

The above map shows that the highest proportion of armed conflicts in the westernmost Rakhine state and easternmost Shan state involves violence against civilians. This state is the epicentre of the 2017 Rohingya crisis facing the eponymous Muslim ethnic group, discrimination against whom, in fact, goes as far back as the 1970s (Albert & Maizland, 2020). In joint second are the regions of Bago and Yangon.

TH1.5.2 Overarching Spatial-Temporal KDE

Next, we can try spatial-temporal kernel density estimation (STKDE) on the entire study area. As always, we start by creating ppp objects out of our sf objects, and checking for duplicates.

acled_mya_sf_quarters <- acled_mya_sf %>% select(Quarter_num)
acled_mya_quarters_ppp <- as.ppp(acled_mya_sf_quarters)
acled_mya_quarters_ppp
Marked planar point pattern: 53653 points
marks are numeric, of storage type  'double'
window: rectangle = [-208804.4, 640934.5] x [1103500.1, 3042960.3] units
any(duplicated(acled_mya_quarters_ppp))
[1] TRUE

Since there are duplicates, we need to handle them using jittering.

acled_mya_quarters_ppp <- rjitter(acled_mya_quarters_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(acled_mya_quarters_ppp))
[1] FALSE

Finally, we join the owin object mapping the borders of Myanmar to the ppp object.

acled_mya_quarters_ppp <- acled_mya_quarters_ppp[mmr_admin1_owin]
plot(acled_mya_quarters_ppp)

We can now perform spatial-temporal KDE (STKDE).

st_kde <- spattemp.density(acled_mya_quarters_ppp)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
summary(st_kde)
Spatiotemporal Kernel Density Estimate

Bandwidths
  h = 42819.86 (spatial)
  lambda = 0.0051 (temporal)

No. of observations
  53559 

Spatial bound
  Type: polygonal
  2D enclosure: [-210008.6, 724647.6] x [1072026, 3158467]

Temporal bound
  [1, 15]

Evaluation
  128 x 128 x 15 trivariate lattice
  Density range: [1.019043e-25, 2.08942e-10]

After performing STKDE, we are able to plot the quarterly STKDE objects between 2021 and September 2024. The full set of plots can be found in Annex A, but the plot for Q1 2021 (id 15) is of particular interest because of the larger than usual densities around the Sagaing and Chin areas, as well as near Yangon.

i = 15
plot(st_kde, quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for", quarters[i]))

TH1.5.3 STKDE on Data Subset

Knowing that a significant part of the armed violence in Myanmar is due to systemic discrimination on the part of the government or military, it would make sense to repeat the process and plot STKDEs on the subset of the data concerning violence against civilians.

[1] "Battles"                    "Explosions/Remote violence"
[3] "Strategic developments"     "Protests"                  
[5] "Violence against civilians" "Riots"                     
acled_sf_againstCivilians <- acled_mya_sf %>% filter(event_type == "Violence against civilians")
# acled_sf_againstCivilians
acled_sf_againstCivilians <- acled_sf_againstCivilians %>% select(Quarter_num)
acled_mya_againstCivilians_ppp <- as.ppp(acled_sf_againstCivilians)
any(duplicated(acled_mya_againstCivilians_ppp))
[1] TRUE
acled_mya_againstCivilians_ppp <- rjitter(acled_mya_againstCivilians_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(acled_mya_againstCivilians_ppp))
[1] FALSE
acled_mya_againstCivilians_ppp <- acled_mya_againstCivilians_ppp[mmr_admin1_owin]
plot(acled_mya_againstCivilians_ppp)

st_kde_againstCivilians <- spattemp.density(acled_mya_quarters_ppp)
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
summary(st_kde_againstCivilians)
Spatiotemporal Kernel Density Estimate

Bandwidths
  h = 42819.86 (spatial)
  lambda = 0.0051 (temporal)

No. of observations
  53559 

Spatial bound
  Type: polygonal
  2D enclosure: [-210008.6, 724647.6] x [1072026, 3158467]

Temporal bound
  [1, 15]

Evaluation
  128 x 128 x 15 trivariate lattice
  Density range: [1.019043e-25, 2.08942e-10]
plot(st_kde_againstCivilians, quarter_numbers[15], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for violence against civilians events in quarter", quarters[i]))

TH1.6 Zooming in

The next step will be to perform spatial point pattern analysis (SPPA) and spatial-temporal point patterns analysis (STPPA). To do this effectively we need to select certain regions to analyse.

While not entirely effective, the preliminary KDE plot identified large clusters of armed conflict in the Yangon area in the south off the Andaman sea, as well as the Sagaing region and Chin state to the central-west. From the chropleth maps, Yangon was shown as the area with the second highest degree of civilian targeting, while Sagaing had the most fatalities from the conflict.

Therefore, it makes sense to focus on Sagaing, Shan, Rakhine and Chin. Additionally, an analysis on the Tanintharyi region will cover the Kra Isthmus, the narrow isthmus linking central Thailand to Peninsula Malaysia, to study the scale of the part of the conflict that is geographically closest to the more stable and prosperous Malaysia and Singapore.

TH1.7 First-order Spatial Point Patterns Analysis

TH1.7.1 Extracting and plotting study areas

The first step is to extract the study areas and create owin objects.

study_names <- c("Sagaing", "Shan", "Rakhine", "Chin", "Yangon", "Tanintharyi")
study_areas <- list("Sagaing" = (mmr_admin1 %>% filter(ST == "Sagaing")),
                    "Shan" = (mmr_admin1 %>% filter(ST == "Shan")),
                    "Rakhine" = (mmr_admin1 %>% filter(ST == "Rakhine")),
                    "Chin" = (mmr_admin1 %>% filter(ST == "Chin")),
                    "Yangon" = (mmr_admin1 %>% filter(ST == "Yangon")),
                    "Tanintharyi" = (mmr_admin1 %>% filter(ST == "Tanintharyi"))
)
par(mfrow = c(1, 4))
for (study_name in study_names) {
  print(paste("Study area:", study_name))
  plot(study_areas[[study_name]], main = study_name) 
}
[1] "Study area: Sagaing"

[1] "Study area: Shan"

[1] "Study area: Rakhine"

[1] "Study area: Chin"

[1] "Study area: Yangon"

[1] "Study area: Tanintharyi"

study_ppps <- list()
study_ppps_original <- list()
study_owins <- list()
par(mfrow = c(1, 3))
for (study_name in study_names) {
  print(paste("Creating ppp for", study_name))
  study_owin <- as.owin(study_areas[[study_name]])
  study_owins[[study_name]] <- study_owin
  study_ppps_original[[study_name]] <- acled_mya_quarters_ppp[study_owin]
  study_ppps[[study_name]] <- rescale.ppp(acled_mya_quarters_ppp[study_owin], 1000, "km")
  plot(study_ppps[[study_name]], main=paste("ppp for", study_name))
}
[1] "Creating ppp for Sagaing"
[1] "Creating ppp for Shan"
[1] "Creating ppp for Rakhine"

[1] "Creating ppp for Chin"
[1] "Creating ppp for Yangon"
[1] "Creating ppp for Tanintharyi"

TH1.7.2 Fixed-bandwidth KDE

fixed_bandwidth_kde_acled <- list()
par(mfrow = c(1, 3))
for (study_name in study_names) {
  fixed_bandwidth_kde_acled[[study_name]] <- density(study_ppps[[study_name]],
               sigma=bw.scott,
               edge=TRUE,
               kernel="gaussian"
               )
  plot(fixed_bandwidth_kde_acled[[study_name]], 
       main = paste("Fixed-bandwidth KDE for", study_name))
}

The fixed-bandwidth KDE returned dense clusters at the following areas for each state or region:

  • Sagaing region: toward the south and southwest
  • Shan state: towards the northwest, with a secondary cluster towards the southwest
  • Rakhine state: one large cluster towards the northwest
  • Chin state: one cluster in the central-north and southeast
  • Yangon region: one dense cluster towards the south coast
  • Tanintharyi region: one dense cluster skirting the northwest coast

TH1.7.3 Adaptive-bandwidth KDE

adaptive_bandwidth_kde_acled <- list()
par(mfrow = c(1, 3))
for (study_name in study_names) {
  adaptive_bandwidth_kde_acled[[study_name]] <- adaptive.density(study_ppps[[study_name]],
               method = "kernel"
               )
  plot(adaptive_bandwidth_kde_acled[[study_name]], 
       main = paste("Adaptive-bandwidth KDE for", study_name))
}

Fixed bandwidth KDE using bw.scott and the gaussian kernel returned a more meaningful result than adaptive bandwidth KDE.

The bw.scott method applies Scott’s Rule for Bandwidth Selection for Kernel Density, which calculates a bandwidth \(\sigma\) that is proportional to \(n^{\frac{-1}{d-4}}\), for number of points \(n\) and number of spatial dimension \(d\). This rule is fast to compute and, as we’ve seen, was able to produce a larger bandwidth than the mean-square error criterion defined with the bw.diggle() method (Diggle, 1985)

Just for completeness, we will now attempt to map the results of fixed-bandwidth KDE for Rakhine state by creating a RasterLayer.

rakhine_kde <- fixed_bandwidth_kde_acled[["Rakhine"]]
gridded <- as(rakhine_kde, "SpatialGridDataFrame")
spplot(gridded)

raster <- raster(rakhine_kde)
raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 2.158983, 3.656088  (x, y)
extent     : -210.0086, 66.34123, 1921.825, 2389.805  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : 0.0007171686, 0.3808087  (min, max)
projection(raster) <- CRS("+init=EPSG:32647")
raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 2.158983, 3.656088  (x, y)
extent     : -210.0086, 66.34123, 1921.825, 2389.805  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=47 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : layer 
values     : 0.0007171686, 0.3808087  (min, max)
tm_shape(raster) + tm_raster("layer", palette = "viridis") + tm_layout(legend.position = c("left", "bottom"), frame = FALSE)

TH1.7.4 Spatial-temporal KDE

We can also plot spatial-temporal KDEs for each of our selected regions. The full plots can be found in Annex B.

stkde_six_regions <- list()
for (study_name in study_names) {
  print(study_name)
  stkde_six_regions[[study_name]] <- spattemp.density(study_ppps[[study_name]])
}
[1] "Sagaing"
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
[1] "Shan"
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
[1] "Rakhine"
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
[1] "Chin"
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
[1] "Yangon"
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.
[1] "Tanintharyi"
Calculating trivariate smooth...Done.
Edge-correcting...Done.
Conditioning on time...Done.

TH1.7.5 Nearest neighbour analysis

For each of the 6 regions we are analysing, we can use nearest neighbour analysis to find out whether the spatial point patterns of our armed conflict events exhibit clustering, dispersion or random distribution. The nearest neighbour analysis is done by performing the Clark-Evans test on the hypothesis that the distribution of armed conflict events is random.

for (study_name in study_names) {
  print(study_name)
  ppp <- study_ppps[[study_name]]
  # plot(ppp)
  print(clarkevans.test(ppp, correction = "none", clipregion = NULL, alternative = c("clustered"), nsim = 99))
}
[1] "Sagaing"

    Clark-Evans test
    No edge correction
    Z-test

data:  ppp
R = 0.16514, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

[1] "Shan"

    Clark-Evans test
    No edge correction
    Z-test

data:  ppp
R = 0.1622, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

[1] "Rakhine"

    Clark-Evans test
    No edge correction
    Z-test

data:  ppp
R = 0.20771, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

[1] "Chin"

    Clark-Evans test
    No edge correction
    Z-test

data:  ppp
R = 0.17494, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

[1] "Yangon"

    Clark-Evans test
    No edge correction
    Z-test

data:  ppp
R = 0.13285, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

[1] "Tanintharyi"

    Clark-Evans test
    No edge correction
    Z-test

data:  ppp
R = 0.12734, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

The nearest neighbour analyses for all our regions returned identical p-values of 2.2e-16. Thus, we can conclude that the spatial point patterns for all our target regions are clustered. We can corroborate this test result with the fixed-bandwidth KDE visualisations, which indicate the presence of densely-packed conflict areas.

TH1.8 Second order spatial point patterns analysis

Second order SPPA investigates variations in observations due to the way they interact with one another. In this section, we will use the F- and G-functions to analyse the spatial point patterns in our 6 selected regions.

Ripley’s K-function valuates the distribution of points over a range of distances, giving insight into spatial clustering or dispersion at different scales. In contrast, the F-function (nearest neighbor distribution) and G-function (empty space function) are limited to analyzing specific distances or single-scale proximity between points. This makes Ripley’s K-function a more comprehensive analysis of spatial point patterns at multiple spatial scales, but supposedly more computationally complex. An attempt to run Ripley’s K-function on armed conflict data in Sagaing - both for the entire study area and just in 2022 - led to a half-hour long freeze, so it will be omitted.

TH1.8.1 Extracting study period

The functions run slower when fed with more data. Hence, we need to extract a smaller study period. Since the year 2022 saw more fatalities, it makes sense to extract only the events that occurred in 2022.

acled_mya_2022_sf <- acled_mya_sf %>% filter(Quarter_num >= 5 & Quarter_num < 7) %>% select(Quarter_num)
acled_mya_2022_ppp <- as.ppp(acled_mya_2022_sf)
acled_mya_2022_ppp <- rjitter(acled_mya_2022_ppp, retry=TRUE, nsim=1, drop=TRUE)
study_ppps_2022 <- list()
for (study_name in study_names) {
  print(paste("Creating ppp for", study_name))
  study_owin <- study_owins[[study_name]]
  study_ppps_2022[[study_name]] <- acled_mya_2022_ppp[study_owin]
}
[1] "Creating ppp for Sagaing"
[1] "Creating ppp for Shan"
[1] "Creating ppp for Rakhine"
[1] "Creating ppp for Chin"
[1] "Creating ppp for Yangon"
[1] "Creating ppp for Tanintharyi"

TH1.8.1 Sagaing Region

TH1.8.1.1 G-function

study_name = "Sagaing"
G_function <- Gest(study_ppps_2022[[study_name]], correction = "border")
plot(G_function, xlim=c(0,500), main = paste("G-function estimation for", study_name))

G_csr <- envelope(ppp, Gest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(G_csr, main = paste("G function test for spatial randomness for", study_name))

TH1.8.1.2 F-function

F_function <- Fest(study_ppps_2022[[study_name]])
plot(F_function, xlim=c(0,500), main = paste("F-function estimation for", study_name))

Error in bestlegendpos(...) : All objects were empty
F_csr <- envelope(ppp, Fest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(F_csr, main = paste("F function test for spatial randomness for", study_name))

TH1.8.2 Shan State

TH1.8.2.1 G-function

study_name = "Shan"
G_function <- Gest(study_ppps_2022[[study_name]], correction = "border")
plot(G_function, xlim=c(0,500), main = paste("G-function estimation for", study_name))

G_csr <- envelope(ppp, Gest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(G_csr, main = paste("G function test for spatial randomness for", study_name))

TH1.8.2.2 F-function

F_function <- Fest(study_ppps_2022[[study_name]])
plot(F_function, xlim=c(0,500), main = paste("F-function estimation for", study_name))

Error in bestlegendpos(...) : All objects were empty
F_csr <- envelope(ppp, Fest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(F_csr, main = paste("F function test for spatial randomness for", study_name))

TH1.8.3 Rakhine State

TH1.8.3.1 G-function

study_name = "Rakhine"
G_function <- Gest(study_ppps_2022[[study_name]], correction = "border")
plot(G_function, xlim=c(0,500), main = paste("G-function estimation for", study_name))

G_csr <- envelope(ppp, Gest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(G_csr, main = paste("G function test for spatial randomness for", study_name))

TH1.8.3.2 F-function

F_function <- Fest(study_ppps_2022[[study_name]])
plot(F_function, xlim=c(0,500), main = paste("F-function estimation for", study_name))

Error in bestlegendpos(...) : All objects were empty
F_csr <- envelope(ppp, Fest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(F_csr, main = paste("F function test for spatial randomness for", study_name))

TH1.8.4 Chin State

TH1.8.4.1 G-function

study_name = "Chin"
G_function <- Gest(study_ppps_2022[[study_name]], correction = "border")
plot(G_function, xlim=c(0,500), main = paste("G-function estimation for", study_name))

G_csr <- envelope(ppp, Gest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(G_csr, main = paste("G function test for spatial randomness for", study_name))

TH1.8.4.2 F-function

F_function <- Fest(study_ppps_2022[[study_name]])
plot(F_function, xlim=c(0,500), main = paste("F-function estimation for", study_name))

F_csr <- envelope(ppp, Fest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(F_csr, main = paste("F function test for spatial randomness for", study_name))

TH1.8.5 Yangon Region

TH1.8.5.1 G-function

study_name = "Yangon"
G_function <- Gest(study_ppps_2022[[study_name]], correction = "border")
plot(G_function, xlim=c(0,500), main = paste("G-function estimation for", study_name))

G_csr <- envelope(ppp, Gest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(G_csr, main = paste("G function test for spatial randomness for", study_name))

TH1.8.5.2 F-function

F_function <- Fest(study_ppps_2022[[study_name]])
plot(F_function, xlim=c(0,500), main = paste("F-function estimation for", study_name))

Error in bestlegendpos(...) : All objects were empty
F_csr <- envelope(ppp, Fest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(F_csr, main = paste("F function test for spatial randomness for", study_name))

TH1.8.6 Tanintharyi Region

TH1.8.6.1 G-function

study_name = "Tanintharyi"
G_function <- Gest(study_ppps_2022[[study_name]], correction = "border")
plot(G_function, xlim=c(0,500), main = paste("G-function estimation for", study_name))

G_csr <- envelope(ppp, Gest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(G_csr, main = paste("G function test for spatial randomness for", study_name))

TH1.8.6.2 F-function

F_function <- Fest(study_ppps_2022[[study_name]])
plot(F_function, xlim=c(0,500), main = paste("F-function estimation for", study_name))

F_csr <- envelope(ppp, Fest, nsim = 99)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(F_csr, main = paste("F function test for spatial randomness for", study_name))

Since all F-functions lie entirely below the envelope, and all G-functions lie entirely above the envelope, we can conclude that the armed conflict areas exhibit clustering. This is backed up by the KDE plots, which show the presence of dense clusters.

TH1.9 Concluding thoughts

In conclusion, we have explored different spatial point pattern analysis methods using armed conflict data in Myanmar. We attempted spatial-temporal KDE on the entire study period and identified dense areas around the Sagaing and Chin areas. Following that, we identified a subset of states and regions for deeper analysis based on exploratory data analysis and the STKDE plots and attempted fixed- and adaptive-bandwidth KDE on these areas, as well as second-order SPPA using the F- and G-functions. We have found bw.scott() to be more effective in performing KDE on this dataset than the usual bw.diggle(), and the conflict points in all regions to exhibit clustering.

Now that we know how the points behave and where they are clustered, one possible extension of this study could be to do deeper research on Myanmar’s municipalities, townships and diverse ethnic groups, and repeat the SPPA process with the admin2 or admin3 plots. This is so that we better identify the needs and desires of each ethnic group, helping diplomats and policymakers formulate more targeted solutions for the ethnic divide and conflict in the country. In addition, Ripley’s K- function, as well as Besag’s L-function, could be attempted again given more time and computing power, along with the incorporation of other data sources, such as economic data, which could enhance our understanding of the relationship between the conflict and the country’s overall economic development. It is worth noting, however, that no amount of geospatial analytics or field studies will be effective without a strong, morally upright and benevolent government in Myanmar, with the political will and know-how to formulate and implement sound policies in response to analyses such as this one, for the betterment of their people and country.

As a long-time Python programmer and learner of data analytics, this exercise has been a rigourous introduction not only to R programming, but to wrangling and analysing points in geographical space. Though frustrating at times, I strongly believe that the best learning comes from experiencing and recovering from adversity, and I look forward to the subsequent exercises and concepts to come.

References

  1. Albert, E., & Maizland, L. (2020, January 23). What forces are fueling Myanmar’s Rohingya crisis? Council on Foreign Relations. https://www.cfr.org/backgrounder/rohingya-crisis
  2. Diggle, P.J. (1985). A kernel method for smoothing point process data. Applied Statistics (Journal of the Royal Statistical Society, Series C), 34, 138–147.
  3. Maizland, L. (2022, January 31). Myanmar’s troubled history: Coups, military rule, and ethnic conflict. Council on Foreign Relations. https://www.cfr.org/backgrounder/myanmar-history-coup-military-rule-ethnic-conflict-rohingya
  4. Myanmar Information Management Unit. (n.d.). Country Overview. Retrieved 16 September 2024, from https://www.themimu.info/country-overview
  5. Raleigh, C., Kishi, R. & Linke, A. Political instability patterns are obscured by conflict dataset scope conditions, sources, and coding choices. Humanit Soc Sci Commun 10, 74 (2023). https://doi.org/10.1057/s41599-023-01559-4
  6. Rieffel, L. (2019, December 6). Peace and war in Myanmar. Brookings. https://www.brookings.edu/articles/peace-and-war-in-myanmar/
  7. Scott, D.W. (1992). Multivariate Density Estimation. Theory, Practice and Visualization. New York: Wiley.
  8. Smith, M. (1999). Burma: Insurgency and the politics of ethnicity (2. (updated) ed). Univ. Press [u.a.].
  9. Sullivan, D. P. (2021, October 21). Dire consequences: Addressing the humanitarian fallout from myanmar’s coup. Refugees International; Refugees International. https://www.refugeesinternational.org/reports-briefs/dire-consequences-addressing-the-humanitarian-fallout-from-myanmars-coup/
  10. UNDP (Ed.). (2024). Breaking the gridlock: Reimagining cooperation in a polarized world. United Nations Development Programme (UNDP).

Annex A: Full Quarterly STKDE Plots for Jan 2021- Sep 2024

A1: Full STKDE Plots

for (i in 1: length(quarter_numbers)) {
  plot(st_kde, quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for", quarters[i]))
}

A2: For events with civilian targeting only

for (i in 1: length(quarter_numbers)) {
  plot(st_kde_againstCivilians, quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("Spatial-temporal KDE for violence against civilians events in quarter", quarters[i]))
}

Annex B: Region-specific STKDE plots

for (study_name in study_names) {
  for (i in 1: length(quarter_numbers)) {
    print(study_name)
    print(quarter_numbers[i])
    plot(stkde_six_regions[[study_name]], quarter_numbers[i], override.par=FALSE, fix.range=TRUE, main=paste("STKDE for", study_name , "in quarter", quarters[i]))
  }
}
[1] "Sagaing"
[1] 15

[1] "Sagaing"
[1] 14

[1] "Sagaing"
[1] 13

[1] "Sagaing"
[1] 12

[1] "Sagaing"
[1] 11

[1] "Sagaing"
[1] 10

[1] "Sagaing"
[1] 9

[1] "Sagaing"
[1] 8

[1] "Sagaing"
[1] 7

[1] "Sagaing"
[1] 6

[1] "Sagaing"
[1] 5

[1] "Sagaing"
[1] 4

[1] "Sagaing"
[1] 3

[1] "Sagaing"
[1] 2

[1] "Sagaing"
[1] 1

[1] "Shan"
[1] 15

[1] "Shan"
[1] 14

[1] "Shan"
[1] 13

[1] "Shan"
[1] 12

[1] "Shan"
[1] 11

[1] "Shan"
[1] 10

[1] "Shan"
[1] 9

[1] "Shan"
[1] 8

[1] "Shan"
[1] 7

[1] "Shan"
[1] 6

[1] "Shan"
[1] 5

[1] "Shan"
[1] 4

[1] "Shan"
[1] 3

[1] "Shan"
[1] 2

[1] "Shan"
[1] 1

[1] "Rakhine"
[1] 15

[1] "Rakhine"
[1] 14

[1] "Rakhine"
[1] 13

[1] "Rakhine"
[1] 12

[1] "Rakhine"
[1] 11

[1] "Rakhine"
[1] 10

[1] "Rakhine"
[1] 9

[1] "Rakhine"
[1] 8

[1] "Rakhine"
[1] 7

[1] "Rakhine"
[1] 6

[1] "Rakhine"
[1] 5

[1] "Rakhine"
[1] 4

[1] "Rakhine"
[1] 3

[1] "Rakhine"
[1] 2

[1] "Rakhine"
[1] 1

[1] "Chin"
[1] 15

[1] "Chin"
[1] 14

[1] "Chin"
[1] 13

[1] "Chin"
[1] 12

[1] "Chin"
[1] 11

[1] "Chin"
[1] 10

[1] "Chin"
[1] 9

[1] "Chin"
[1] 8

[1] "Chin"
[1] 7

[1] "Chin"
[1] 6

[1] "Chin"
[1] 5

[1] "Chin"
[1] 4

[1] "Chin"
[1] 3

[1] "Chin"
[1] 2

[1] "Chin"
[1] 1

[1] "Yangon"
[1] 15

[1] "Yangon"
[1] 14

[1] "Yangon"
[1] 13

[1] "Yangon"
[1] 12

[1] "Yangon"
[1] 11

[1] "Yangon"
[1] 10

[1] "Yangon"
[1] 9

[1] "Yangon"
[1] 8

[1] "Yangon"
[1] 7

[1] "Yangon"
[1] 6

[1] "Yangon"
[1] 5

[1] "Yangon"
[1] 4

[1] "Yangon"
[1] 3

[1] "Yangon"
[1] 2

[1] "Yangon"
[1] 1

[1] "Tanintharyi"
[1] 15

[1] "Tanintharyi"
[1] 14

[1] "Tanintharyi"
[1] 13

[1] "Tanintharyi"
[1] 12

[1] "Tanintharyi"
[1] 11

[1] "Tanintharyi"
[1] 10

[1] "Tanintharyi"
[1] 9

[1] "Tanintharyi"
[1] 8

[1] "Tanintharyi"
[1] 7

[1] "Tanintharyi"
[1] 6

[1] "Tanintharyi"
[1] 5

[1] "Tanintharyi"
[1] 4

[1] "Tanintharyi"
[1] 3

[1] "Tanintharyi"
[1] 2

[1] "Tanintharyi"
[1] 1